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I.  INTRODUCTION 


The  study  of  direct  numerical  simulation  of  steady  viscous  flow  due  to  a  sphere 
which  rotates  about  a  diameter  with  constant  angular  velocity  in  a  fluid  at  rest  has 
historically  been  the  focus  of  a  great  deal  of  scientific  attention.  Most  of  the  previous 
studies  were  based  on  finite-difference  schemes. 

In  this  study,  a  different  numerical  scheme  based  on  spectral  methods  was  used. 
The  validity  of  the  numerical  method  and  the  results  were  checked  with  the  previous 
studies.  Spectral  methods  have  become  increasingly  popular  in  recent  years,  especially 
since  the  development  of  fast  transform  methods,  with  applications  in  numerical  weather 
prediction,  numerical  simulation  of  laminar  and  turbulent  flows  and  other  problems  where 
high  accuracy  is  desired  for  complicated  results. 

The  objective  of  this  study  is  to  compute  the  numerical  simulation  of  a  steady  flow 
about  a  spinning  sphere  by  using  the  spectral  numerical  methods.  A  Matlab  program  was 
developed  by  using  spectral  methods.  However  the  method  itself  was  not  compared  to  the 
other  methods  used  in  previous  studies.  The  reliability  and  validity  of  the  numerical 
method  were  investigated  to  improve  our  current  ability  to  compute  and  predict  the 
evolution  of  flow  for  a  particular  geometry  and  conditions. 


II.  BACKGROUND  STUDIES 


The  spectral  method  is  based  on  the  collocation  method  which  appears  to  have 
been  used  first  by  Slater  and  Kantorovic  in  specific  applications  in  1934.  It  was  developed 
as  a  general  method  for  solving  ordinary  differential  equations.  This  method  was  revived 
by  Wrigt  in  1964.  The  applications  of  Chebyshev  polynomial  expansions  to  the  initial 
value  problems  were  involved  in  these  studies. 

The  spectral  collocation  method  was  applied  to  partial  differential  equations  for 
spatially  periodic  problems  by  Kreiss  and  Oliger,  who  called  it  the  Fourier  method,  and 
Orszag,  who  termed  it  "pseudospectral".  These  were  the  earliest  applications  of  spectral 
collocation  or  the  pseudospectral  method.  This  approach  was  very  attractive  because  of 
its  application  to  variable-coefficient  and  even  non-linear  problems. 

The  Galerkin  approach  which  depends  on  the  same  trial  and  test  functions,  was 
applied  to  a  meteorological  model  by  Silberman  in  1954.  This  was  the  first  serious 
application  of  spectral  methods  to  PDEs.  In  1970,  Machenhauer  and  Rasmussen 
developed  transform  methods  for  evaluating  convolution  sums  arising  from  quadratic  non- 
linearity.  Spectral  Galerkin  methods  became  practical  for  high  resolution  calculations  of 
such  non-linear  problems.  Applications  in  fluid  dynamics  were  reviewed  in  the  symposium 
proceedings  edited  by  Voigt,  Gottlieb  and  Hussaini  in  1984. 

For  the  problem  of  the  flow  induced  by  a  spinning  sphere  being  considered  in  this 
study,  the  flow  behavior  for  small  values  of  the  Reynolds  number,  Stokes  [Ref.  1],  Lamb 
[Ref.  2],  Bickley  [Ref  .3],  Collins  [Ref.  4],  Thomas  and  Walters  [Ref.  5],  Ovseenko 
[Ref.  6],  and  Takagi  [Ref.  7]  have  given  theoretical  results.  For  high  Reynolds  number,  the 
boundary-layer  equations  have  been  studied  by  Howarth  [Ref.  8],  Nigam  [Ref.  9], 
Stewartson  [Ref.  10],  Fox  [Ref.  11],  Banks  [Refs.  12  &  13],  Manohar  [Ref.  14],  and 
Singh  [Ref.  15].  Over  a  large  range  of  values  of  Reynolds  number  the  flow  has  been 
investigated  experimentally  by  Kobashi  [Ref.  16],  Bowden  and  Lord  [Ref.  17],  Kreith  et 
al  [Ref.  18],  and  Sawatzki  [Ref.  19].  Another  numerical  method  based  on  the  use  of 


specialized  techniques  to  obtain  an  approximation  of  the  finite-difference  form  of  the  full 
partial  differential  equations  was  performed  by  Allen  and  Southwell  [Ref.  20],  Dennis 
[Refs.  21  &  22],  Roscoe  [Refs.  23  &  24],  Spalding  [Ref.  25]  and  Dennis,  Ingham  and 
Singh  [Ref.  26]. 


A. 


III.  GOVERNING  EQUATIONS 

DERIVATION  OF  EQUATIONS 

The  governing  equations  for  the  problem  are  as  follow; 


3u         ^  1  „        „, 

— +  (u.V)u  =  --Vp  +  vV2u 
dt  p 


V.u  =  0 


Define  dimensionless  variables; 


r    =—  u    = —  Vp    = - 

a  U  y       pU2 


where     a  :  radius  of  the  sphere, 
U  :  free  stream  velocity, 
p   :  density  of  the  fluid, 
v   :  kinematic  viscosity  of  the  fluid. 

Introduce  the  dimensionless  variables  into  Equation  (3.1); 


(3.1)  Momentum  equation 


(3.2)  Continuity  equation 


a 


Uv  du*      U2      .  „     . 
—— +  — (u  .V)u  = 
a"    dt        a 


1    U2T7     *  ^ 

Vp  +v-rV2u 

pa  a 


Divide  Equation  (3.3)  by 


U.v 


du*      U.a      .  „     , 

—  + (u  .V)u 

dt         v 


U.a 


v 


Vp'+vVV 


(3.3) 


(3.4) 


Hereafter  the  superscript  *  will  be  omitted  for  nondimensional  terms.    Define  the 

U(2a) 

Reynolds  number  based  on  the  diameter,  Red  =  ,  and  take  the  Curl  of  Equation 

v 

(3.4)  to  be  able  to  eliminate  the  pressure  term  and  introduce  vorticity,  co  =  V  x  u    where 

5 


CO=irOV+ieC0e+id>CL)kt>. 


The  transport  equation  for  the  O-component  of  the  vorticity  is 


—  gv0-Vx(ux co 0i0)  =  - —  V--    ,    .   , 


at 


Re 


r"  sin"  By 


CO^lo 


(3.5) 


Convert  the  velocity  in  Equation  (3.5)  by  defining  the  axisymmetric  Stream  function,  V|/, 
and  an  "angular  circulation",  Q. 


u  =  Vx 


M> 


rsin0  ° 


Q 


+  1 


*  rsin6 


<j)  ;  Azimuthal  direction     (3.6) 


In  spherical  coordinates,  Equation  (3.6)  is  obtained  as 


u  = 


1        d\\f 

r2sin6  36 


l    - 


d\\f 


rsinw  or 


6  di 


,e+—- "I, 
rsin9 


(3.7) 


where 


u  =  urir +ueie+u«,i«) 


(3.8) 


In  the  free  stream  \it  =  —  Ur2  sin2  0  or  \i/=  —  r2sin20  in  nondimensional  form. 

Y      2  Y      2 

By  using  the  definition  of  the  stream  function  it  can  be  shown  that 


Git>  = 


a2\|/     i    i    a  (  i   d\\f 


+ 


dr2       r2  sinG  36 


sine  ae 


rsinG 


(3.9) 


Define  operator  D"  as 


— ,       d2      sin0  3 


dr2  '    r2    90 


(    1      d\ 

sineaey 


(3.10) 


to  obtain  the  azimuthal  component  of  the  vorticity; 


co*  =  - — — ~D2\\f 
rsin6 


(3.11) 


Call  [  V  x  (u   x  con)i<j))].i(t)    in  Equation  (3.5)  nonlinear  Term  G(r,9,t) .  By  carrying  out  the 
calculations  in  spherical  coordinates 


1 


Vx(uxo),iJ=  — 

r   sin0 


3(\i/,D2\i/)       _, 

ibr+2D¥LW 


(3.12) 


is  obtained,  where 


(I  =  COS0 

L  = 


\l     a     1  a 

+— 


l-|i2  3r     r  d\i 
a(z,,z2)      dzt  dz2      3z,  9z2 
a(x,y)        3x    dy       dy    dx 


(3.13) 


Combine  all  terms  to  get 


a   -2         Red    1 


a(\i/,D2\i/)     — , 


=  D> 


(3.14) 


Equation  (3.14)  is  called  the  Vorticity-Stream  function  form  of  the  Momentum 
Equation  in  axisymmetric  spherical  coordinate  system.  For  computational  purposes, 
which  will  become  clear  later,  define  a  modified  stream  function  C(r,6,t),  which  is  related 
to  the  usual  Stokes  stream  function  (\j/)  by  the  relation; 

V|/  =  rCsin0  (3.15) 

Also  for  any  function  f  (r,0) ,  define  operator  D2  such  that 

D2(frsin0)  =  rsin0D2f 


where 


D2=V2- 


r   sin"  0 


This  gives 


03£b=-D2C  (3.16) 

v       rsin0 


u  =  Vx(CiJ+— —  i*  (3.17) 


1       ^ 
rsin0  d0 

1  a 

where  ue=--— (rC)  (3.18) 


r  dr 
Q 
U°_rsin0 

The  modified  stream  function  C  is  written  as  the  sum  C  =  C  +  c ,  where  C  =  0 , 
(which  is  the  prescribed  quiescent  free  stream  condition),  and  c  represents  the  disturbance 
produced  by  the  presence  of  the  spinning  sphere. 


Returning  to  Equation  (3.5),  rearrange  by  multiplying  both  sides  with  r 


2Red 


ReH  3co       ,  Re 


2   at 


2      wd 


=  r 


2    iw,d 


G(r,e,t)  +  r2D"co 


r2D2C  =  -r2co 


2    ' 
(3.19) 
(3.20) 


where 


G(r,e,t)  =  [Vx(uXGV*)]-io 


(3.21) 


Note  that  the  operator  D2  is  defined  as 


D2  = 


r2  3r 


r ,  ^ 


i 


+  , . , 


3r  J     r~  sin2  0  36 


d  f  3^ 

sine  — 

ae 


V 


1 


1       .      ~> 


r"  sin"  6 


(3.22) 


Divide  r2D2  in  two  parts  which  are  D2  and  D2  respectively  and  obtain 


r2D2=D2  + 


D; 


r  '  sin2  6     e 
where  D2  and  D^are  as  defined  in  Chapter  IV. 

The  nonlinear  term  in  Equation  (3.21),  G(r,6,t),  consists  of  8  terms.  These  are 


(3.23) 


i  aco^  ac 
r  3r  ae 


i  aco0  ac 
7  ae  ar 


cote     ac 

r  dr 


CotO    aw, 


ar 


cop  ac 
r2  ae 


c  aco^ 
r2  ae 


2u. 


cote  au^ 

r       ar 


2u$  au0 
r2    ae 


Time  integration  of  vorticity  equation  (3.19)  is  accomplished  through  the  use  of  an 
explicit  second-order  Adams-Bashforth  scheme  for  the  nonlinear  terms  and  an  implicit 


second-order  Crank-Nicolson  scheme  for  the  viscous  and  linear  terms.  The  calculations 
are  made  in  physical  r-space  and  spectral  6-space.  If  we  denote  vorticity  oo(r,t)  as  the 
vector  of  N  sine  series  then  the  discretized  form  of  Equation  (3.19)  is 


Re, 


2    ixv-d 


«t+A,  -©, 


At 


Re 


=  r^[aGI+PG^]+[D;+A] 


CO.  +  C0 


t+At 


L i        i-  —  \-at  j       l—  f     "    "  "Jl  O 

where  a  and  [3  are  suitably  chosen  weighting  parameters  in  the  implicit  scheme. 
Rearrange  Equation  (3.24)  and  obtain 


(3-24) 


At 


CO 


l+At 


At 


CO 


t-r2Red[cxGt+pGt_Al]      (3.25) 


where  a=—,fi  =  —    and    G  =  (V  x  (u  x  C0oio ))  ■  i0 

[d2  +  AJc(r,t  +  At)  =  -r2co(r,t  +  At) 


(3.26) 


We  find  the  governing  equation  for  u<p  in  same  way: 


dQ.        1 
dt       r" 


aOF.o) 


a(r,cose) 


D2a 


Re 


where     r^u,.  and    C  = - 

rSin0        °  rSin9 


(3.27) 


and  introduce    D'Q  =  D~(uorSin0)  =  rSin6D2u0  and  divide  each  side  by  (r  sinO) 
If  Equation  (3.27)  is  rearranged,  the  result  is 


at 


d(u9,C) 
9(r.6) 


+  CM 


V  C 


^,Ln(rSin6) 


M 


d(r,9) 


+  -^D"u0       (3.28) 
ReH 


10 


The  nonlinear  term  in  Equation  (3.28)  is  called  H(r,0,t)  and  consists  of 


1    du0    dC 
r     dr      30 


1    9u«,    dC 
r     38      dr 


cote     ac 

~u*^r~ 


cote    du^ 
c — — 

r  dr 


r2  ae 


C    du 


<\> 


2   ae 


up  is  written  as  the  sum  of  two  parts:  potential  and  disturbance  parts, 


_-  -        Sin0 

Up  —  Up  +11^         where         u<t  = — r-    is  the  Stokes  flow 


solution  which  satisfies   D"  Up  =  0 .  The  final  form  of  Equation  (3.28)  is 


*i_H(,iM+JLD.W 


(3.29) 


u,(t  +  AO-u,(t)^[3H(t)_H(t_At)3  +  ^ 


UpW  +  u^t  +  At) 


(3.30) 


u«,(t  +  At) 


r  Re,      r2D 


2r^2 


2At  2 


=  u0(t) 


r2Re,      r2D2 


■+■ 


2At  2 


r2Re, 


+  ■ 


[3H(t)-H(t-At)](3.31) 


Re.r2 

Dr2+A-- 


At 


u0(t  +  At)  =  - 


r2  Red 

Dr2+A+- 


At 


r2  Rp 

u0(t)-^-p-[3H(t)-H(t-At)](3.32) 
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B.         BOUNDARY  CONDITIONS 

The  boundary  conditions  for  the  numerical  scheme  are 

3c        3C 
c  =  -C  =  0         —  =  -—  =  0  uo=0      at  r=l  (3.33) 

dr         dr 

c  =  0       co=0     ^  =  0    atr  =  r„  (3.34) 

The  Neumann  Boundary  Conditions  in  Equation  (3.33)  are  handled  with  the 
Green's  Function  Method  that  will  be  explained  later.  Zero  boundary  conditions  in  the  0 
direction  are  automatically  satisfied  with  the  choice  of  the  sine  expansions. 


12 


IV.    METHOD  OF  SOLUTION 


A.         NUMERICAL  METHOD 

A  pseudospectral  or  collocation  method  based  on  the  technique  developed  by 
Marcus  &  Tuckerman  [Ref.  29],  is  used  to  solve  the  equations  of  motion  in  time.  Spectral 
methods  may  be  viewed  as  an  extreme  development  of  weighted  residual  (MWR)  which  is 
generally  used  in  a  discretization  scheme  for  differential  methods.  The  trial  functions  and 
test  functions  are  chosen  as  infinitely  differentiable  global  functions.  This  is  one  of  the 
features  which  distinguish  spectral  methods  from  finite-element  and  finite-difference 
methods.  A  pseudospectral  transform  method  is  used  in  this  study.  The  approach  taken  in 
the  transform  method  is  to  use  the  Fast  Fourier  Transform  (FFT)  to  transform  the 
functions  to  spectral  space  or  Inverse  Fast  Fourier  Transform  (IFFT)  to  transform  the 
equations  to  physical  space. 

A  physical  process  can  be  described  either  in  time  domain,  by  the  values  of  the 
same  h  as  a  function  of  t,  e.g.,  h(t),  or  else  in  frequency  domain,  where  the  process  is 
specified  by  giving  its  amplitude,  H,  as  a  function  of  frequency  f,  that  is  H(f),  with  -«>  <  f  < 
<>=>.  So  once  it  is  known  whether  the  function  is  either  odd  or  even,  this  information  can  be 
used  in  transformations.  If  h(t)  is  real  and  even  then  H(f)  is  real  and  even,  if  h(t)  is  real 
and  odd  then  H(f)  is  imaginary  and  odd.  It  is  necessary  to  go  back  and  forth  between 
these  two  domains  by  means  of  Fourier  transform  equations.  As  explained  below,  odd 
functions,  using  sine  series  expansions  for  CO,  c  and  u<u,  are  used.  The  same  properties  are 
valid  for  derivatives  of  these  expansions.  Using  these  properties  will  increase 
computational  efficiency. 

Functions  are  presented  both  in  spectral  space  as  a  finite  series  of  basis  functions 
and  by  values  at  collocation  grid  points  in  physical  space.  The  vorticity,  CO,  the  modified 
stream  function,  C,  and  the  velocity  in  the  O  direction,  u0,  are  expanded  as  Chebyshev 
polynomials  in  the  radial  direction  and  as  sine  series  in  the  transverse  direction.  The 
products    and   the    derivatives  in  radial  direction  are  done  in  physical  space  while  the 
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derivatives  in  transverse  direction  are  obtained  in  spectral  space.    Since  each  sine  term  in 
the  expansion  satisfies  the  homogenous  9  boundary  conditions  (co  =  0,  4*  =  0,  u<d  =0  at 
9=0,7c)  exactly,  no  further  6  boundary  conditions  need  to  be  applied. 
In  general  let  a  variable,  say  co,  be  expressed  as 

N 

O)(r,0,t)  =  Xa)n(Z't)Sin(0n)  (4.1a) 

n=i 

with 

M 

co„(z,t)  =  £comn(t)Tm(z)  (4.1b) 

m=0 

where  Tm(z)  is  Chebyshev  Polynomial  with  - 1  <  z  <  1 

e.=^7T  "-1A-..N  (4-2) 

z  =  Cos m  =  0,l,...M  (4.3) 

v  m  ; 

An  algebraic  map  is  used  to  map  the  radial  interval  1  <  r  <  r„  to  the  interval 

- 1  <  z  <  1  where 

1  +  z 

r  =  l  +  L ,    lzl<l  (4.4) 

b-z 

2L 

with        b  =  1  + (Note  that  r0  =  1 ,  which  is  the  surface  of  the  sphere.) 

r    -  r„ 

L  is  a  scaling  parameter  which  is  used  to  map  the  radial  interval  to  the  z-interval. 
A  finite,  but  large  outer  r«„,  was  chosen  to  avoid  regularity  problems  in  the  radial 
differentiation. 

The  sine  expansions  (4.1a)  for  co,  C  and  u<t>  satisfy  the  homogenous  G  boundary 
conditions  and  match  exactly  the  periodicity  and  symmetry  conditions  in  6.  The  solution 
of  the  elliptic  equation  (3.20)  is  required  to  obtain  the  modified  stream  function,  C,  from 
vorticity  coat  any  time  level.  Following  [Ref.  29],  separable  derivative  operators,  Dr2  and 
De2,  are  defined  as 
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where 


r2D2=D2+  — ^D 
sin"  0 


df  ,   d\ 


(4.5) 


Dz  =  — 
r      dr^ 

D2  = 


n  a  r 

sin0  — - 


■"SeJ-1 


(4.6) 


In  spectral  0-space  the  effect  of  operator  De2  at  any  fixed  radial  location  (physical  r-space) 
can  be  written  as 


D2.f(0)  = 


'        a2  a 

sin ~  0  — — r  +  sin  0. cos 0 -r—  -  1 

a©2  a© 


IfjSin(je) 


...=  If, 


j=i 


-(l  +  -j2)sin(j0)  +  -j(l  +  j)sin(20  +  j0)  +  -j(l-j)sin(20-j0) 


(4.7) 


S„  = 


(4.8) 


If  this  expression  is  written  in  the  form  of  Dgf(9)  =  \Sl}  fj 

i  =  j  +  2 

i  =  j 

i  =  j-2 

otherwise 

is  obtained. 

This  matrix  is  NxN  if  the  representation  for  D%f  is  truncated  to  N  Fourier  terms. 

Similarly,  the  operation  of  multiplying  a  function  f(9)   by  s'm~2(6)  produces  another 
NxN  matrix  which  is  called  A  matrix  where  the  only  non-zero  terms  are 


f(i-2)(i- 

1) 

4 

l 

<             2 

(i  +  2)(i  + 

1) 

4 

0 

'al(i)  =  -i(i  +  l) 


Aij  "  Ia2(i)  =  -2i 


i  <  j,i  +  j  =  even 


(4.9) 
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Time  integration  of  the  equations  are  obtained  by  using  the  scheme  which  is 
explained  in  Chapter  III. 

The  nonlinear  terms,  G  and  H,  have  to  be  carefully  handled  by  transforming  the 
spectral  coefficients  to  physical  space  first,  then  multiplying  the  terms  in  physical  space, 
and  finally  transforming  the  result  of  these  products  back  to  spectral  space. 

The  radial  derivatives  are  evaluated  by  collocation  methods  and  expressed  as 
matrix  operations  on  the  vector  of  function  values  at  the  grid  points  in  physical  r-space. 
For  (M+l)  collocation  points,  we  introduce  the  (M+l)x(M+l)  matrix  operation,  that  is 


D(l)(n)  =  D;  + 


Re      ' 
a(.)--r' 


(4.10) 


2      2/ 


where,  I  is  the  identity  matrix.   The  operator  [Dr  -r"(Re/At)+A]  in  (3.25)  and  (3.32)  may 
now  be  written  in  block  matrix  form  as 


D(l)(l)          0  a2(l)I  0 

0  D(,)(2)         0  a2(2)I 

0              0  D(1)(3)  0 

0              0              0  D(1)(4) 


(4.11) 


This  matrix  acts  on  the  vector  of  a  length  of  (M+l)xN  for  vorticity  coefficients. 
Equations  (3.25)  and  (3.26)  are  upper  triangular,  block  matrix  problems  in  physical  r- 
space  and  spectral  6-space  that  can  be  solved  by  any  solver  with  the  Dirichlet  and  the 
Neumann  boundary  conditions  discussed  below.  The  equation  (3.32)  is  treated  in  the 
same  manner  and  leads  to  an  upper  triangular,  block  matrix  problem. 
B.        GREEN'S  FUNCTION  METHOD 

To  enforce  the  radial  boundary  conditions  in  c-co  equations,  (3.25)  and  (3.26), 
Green's    Function    Method    has    been    developed    which    gives    the    correct    boundary 
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conditions,  (3.33)  and  (3.34),  and  allows  the  surface  vorticity  to  develop  naturally.    If  a) 
and  c  consist  of  two  parts,  homogenous  and  particular  solutions,  the  following  is  written: 

N  N 

CO  =  COp+X^j«j  and  C  =  CP+X^JCJ  j=l,2,...,N  (4.12) 

j=i  j=i 

If  (4.12)  is  introduced  into  (3.25)  and  (3.26),  the  homogenous  parts  as  follow  with 
corresponding  boundary  conditions  occur. 

E2coj(r,G)  =  0  (4.13)  H2cj=-r265j  (4.14) 

coj(r  =  l,91)  =  8ij  cj(r  =  l,ei)  =  0 

(Oj(r  =  oo,0.)  =  o  cJ(r  =  oc,e,)  =  0 

The  particular  parts  with  corresponding  boundary  conditions  are 


E2cop(r,9,t  +  At)  =  R(r,e,t)   (4.15)        H2cp(r,6,t  +  At)  =  -r2G)p(r,6,t)  (4.16) 


oop(r  =  l,0)  =  O  cp(r  =  l,6)  =  -Clr= 

cop(r  =  oo,9)  =  0  Cp(r  =  <-,6)  =  0 


E2  and  H2  are  obtained  from  the  block  matrix  defined  in  (4.11) 
To  find  A, 


3c        ac 

9rlr='-~ar 


— U=-T-U  (4-17) 


is  enforced. 
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If  (4.12)  is  substituted  into  (4.17), 


'r=i.e,  "r Lu  ;\r  'r=i.e,  ^j 


3r 


ti  dr 


dC 

dr 


'r=l. 


(4.18) 


is  obtained. 


Allow 


3C      3c, 


dv       dr 


-lr  =  ).6, 


5eJ 

=  H,         and         ^    lr=ie,  -Ag 


(4.19) 


Finally  from  (4.18)  and  (4.19),  Hj  =  /.A^X^  occurs  and  we  solve  for  Xj  as  X  =  A"'H. 
These  X  values  can  now  be  used  in  Equation  (4. 12)  to  arrive  at  the  values  of  co  and  c. 
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V.    RESULTS  AND  DISCUSSION 


Numerical  calculations  were  performed  for  following  parameters. 


Re 

M 

N 

L 

r„ 

dt 

2 

64 

64 

7 

150 

0.05 

4 

64 

64 

6 

140 

0.05 

10 

64 

64 

4 

130 

0.03 

20 

64 

64 

4 

125 

0.05 

40 

64 

64 

4 

125 

0.05 

60 

64 

64 

4 

100 

0.03 

80 

64 

64 

3 

70 

0.03 

100 

64 

64 

2 

40 

0.01 

150 

64 

64 

2 

30 

0.03 

175 

64 

64 

2 

25 

0.03 

200 

64 

64 

2 

20 

0.05 

1000 

64 

64 

1 

17 

0.01 

2000 

64 

64 

0.6 

14 

0.01 

3000 

120 

128 

1.5 

15 

0.01 

5000 

128 

128 

1 

8 

0.01 

Table  1.  Parameters  Used  in  Numerical  Calculations 

One  of  the  main  problems  associated  with  many  investigations  that  deal  with 
complete  numerical  solution  of  Navier  Stokes  equations  is  the  behavior  of  numerical 
scheme  at  high  Reynolds  numbers.  It  was  found  here  that  numerical  results  could  be 
obtained  for  higher  values  of  Reynolds  numbers.  For  computational  purposes,  the  results 
for  Reynolds  numbers  up  to  5000  are  presented. 
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For  comparison,  T  (Torque)  is  defined  as 

T=  Jo^r'SinedS* 


(5.1) 


S.r=l 


If  (5.1)  is  nondimensionalized,  the  following  occurs: 


_^-  =  Jc^Sin2ede       ,  where      ar,=r- 


3u0     u^ 
3r   "    r 


'r=l 


(5.2) 


The  small  mapping  parameter,  L,  with  increasing  Reynolds  number  provides  a  finer 
grid  in  the  boundary  layer.  The  re-circulation  of  the  flow  near  the  equator  is  another 
important  result  of  a  spinning  sphere  which  will  be  explained  later.  Since  most  of  the  sign 
changes  in  velocity  profiles  occur  in  this  region,  decreasing  L  much  below  one  is  not  a 
solution  to  gain  some  additional  grid  points  in  the  boundary  layer  for  high  values  of 
Reynolds  number  because  finer  grid  near  the  sphere  means  coarser  grid  in  the  outer 
region. 

Several  solutions  for  large  Reynolds  numbers  have  been  obtained  by  integrating 
numerically  the  boundary  layer  equations  found  in  [Ref.  13]  which  gives 


6.48  a2co 

^77,  where   Rea  =  — 


M  =  — — f77,  where   Rea  =— —  (5.3) 


There  are  some  other  investigations  that  recommend  a  different  numerical  factor 
such  as  5.95  [Ref.  8],  6.54  [Ref.  12]  and  6.53  [Ref.  14].  In  [Ref.  26],  Equation  (5.3)  was 
modified  as 


6.48        31 


The  Constants  6.48  and  31  in  Equation  (5.4)  are  for  Reynolds  numbers  based  on  radius; 
these  constants  should  be  modified  to  9.16  and  62  respectively  for  Reynolds  numbers 
based  on  diameter  as  in  this  study. 
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Red 

Takagi 
[Ref.  7] 

Dennis, 

and  Singh 

[Ref.  20] 

Banks 
[Ref.  13] 

Equation 

(5.4) 

Dennis,  and 

Singh 

[Ref.  26] 

Present 

2 

50.307 

50.309 

6.48 

37.480 

50.305 

50.31 

4 

25.216 

25.218 

4.58 

20.082 

25.216 

25.23 

20 

5.401 

5.399 

2.05 

5.149 

5.398 

5.40 

40 

3.048 

1.45 

2.999 

3.048 

3.02 

100 

1.554 

0.916 

1.496 

1.554 

1.54 

200 

0.966 

0.648 

0.958 

0.966 

0.96 

1000 

0.290 

0.352 

0.348 

0.35 

2000 

0.205 

0.236 

0.234 

0.23 

3000 

0.167 

0.188 

0.19 

5000 

0.129 

0.142 

0.14 

Table  2.  Nondimensional  Torque  (M)  at  Different  Reynolds  Numbers 
Figure  1  shows  the  variation  of  nondimensional  torque  with  Reynolds  numbers.   It 
is  seen  that  as  the  Reynolds  number  increases  the  general  trend  is  quite  consistent  with  the 
boundary  layer  solution  given  in  [Ref.  13]. 


Figure  1.  — ,  Banks  [Ref.  13];  o,  Sawatzki  [Ref.  19];  +,  Present 
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Figure  2  and  Figure  3  show  the  transverse  and  azimuthal  components  of  the 
dimensionless  skin  friction  on  the  sphere,  respectively.  As  it  is  in  Figure  1,  the  higher 
Reynolds  number,  the  closer  it  is  to  the  boundary  layer  solution.  Transverse  and  azimuthal 
components  of  dimensionless  skin  friction  are  defined  as 


1/2 


T0  = 


Re 


90 


dj 


dv 


A=i 


(5.5) 


**  = 


VYao 


VRed  J 


I  5r  A= 


(5.6) 


Figure  2.  Transverse  Component  of  Skin  Friction  ;  — ,  Banks  [Ref.  13] 
o,  Red=20;  +,  Red=200;  *,  Red=3000 
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Figure  3.  Azimuthal  Component  of  Skin  Friction;  — ,  Banks  [Ref.  13] 
o,  Red=20;  +,  Red=200;  *,  Red=3000 

Present  results  agree  well  with  the  boundary  layer  solution  near  the  pole  but  this  is 
not  so  for  the  equator  region.  This  is  because  the  boundary  layer  solution  is  obtained  by 
integrating  parabolic  partial  differential  equations  from  9=0  with  known  initial  conditions 
and  the  solution  near  the  equator  is  determined  from  this  procedure.  Thus,  it  does  not 
satisfy  the  physical  boundary  conditions. 

Since  a  finite  r^  is  used  and  the  modified  stream  function  is  forced  to  be  zero  at  r«, 
and  on  the  sphere,  a  spurious  re-circulating  region  of  extremely  weak  flow  (relative 
magnitude  ~  10~7)  is  introduced  to  account  for  the  continuity  in  the  flow  domain.  The 
inflow  at  the  pole  changes  to  an  outflow  at  the  equator.  The  angular  position  of  this 
change  depends  on  Re  but  does  not  vary  greatly  in  terms  of  radial  distance  from  the 
sphere.  However,  the  general  trend  is  to  get  closer  to  the  sphere  with  increasing  Reynolds 
numbers.  Figure  4  shows  the  stream  lines  in  symmetric  plane  for  different  values  of 
Reynolds  numbers. 
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Figure  4.  Stream  Lines  for  Various  Re 

In  Figure  5,  the  jet  effect  can  be  easily  seen  at  the  equator  for  two  different 
Reynolds  numbers.  The  darker  color  represents  the  regions  where  the  magnitude  of 
dimensionless  velocity  in  the  symmetric  plane  is  higher.  The  reference  velocity  in  Figure  5 
is  the  velocity  in  the  azimuthal  direction  on  the  sphere  which  has  a  value  of  unity. 


0.2 

0.1 

5 

0.1 


Re=200 


Figure  5.  Jet  Effect  at  Equator 
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Figure  6a  and  Figure  6b  show  the  experimental  flow  visualization  obtained  by 
F.  P.  Bowden  and  R.  G.  Lord  [Ref.  17],  and  present  numerical  solution  for  Red=3300 
respectively.  Note  that,  unlike  Figure5,  the  velocity  fields  in  Figure  6a  and  Figure  6b  are  in 
3-D.  As  it  is  seen  there  is  a  reasonable  resemblance  between  experimental  and  numerically 
obtained  pictures.  In  the  experiment,  the  sphere  was  levitated  by  means  of  a  magnetic 
field.  It  was  noted  in  [Ref.  17]  that,  because  of  this  magnetic  field,  there  is  a  slight  upward 
convection  current  so  the  smoke  rises  slowly;  this  explains  why  the  top  half  of  the  picture 
is  white  and  the  bottom  half  black.  Note  also  that,  the  length  of  the  radial  jet  is  relatively 
shorter  in  Figure  6a.  This  behavior  can  be  explained  by  the  presence  of  possible  shear 
turbulence  in  the  radial  jet  for  high  Reynolds  numbers  causing  the  laminar  jet  behavior  to 
disintegrate.  Since  the  present  results  are  obtained  by  means  of  a  numerical  solution 
without  any  kind  of  turbulence  modeling,  there  is  no  effect  of  turbulence  in  our  cases, 
even  for  higher  Reynolds  numbers. 


Figure  6a.  Smoke  Picture  of  the  Boundary  Layer  Flow  for  Re=3300.  (From  Ref.  17) 


Figure  6b.  Numerical  Picture  of  the  Boundary  Layer  Flow  for  Re=33O0. 
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In  the  boundary  layer  approximation,  it  is  assumed  that  the  velocity  profile  in  the 
boundary  layer  just  before  the  equator  is  the  same  as  the  one  at  the  start  of  boundary 
collision.  This  assumption  yields  zero  radial  velocity  at  6  =  n/2,  near  the  sphere.  From 
the  results  it  is  seen  that  radial  velocity  is  not  zero  at  the  equator  near  the  sphere.  Figure  7 
shows  the  variation  of  radial  velocity  with  9  for  Re  =  200  at  different  radial  distances  from 
the  sphere.  Note  that  the  maximum  value  of  radial  velocity  occurs  exactly  at  the  equator. 
Besides,  there  is  an  inward  flow  at  both  poles.  Although  the  maximum  value  of  radial 
velocity  at  the  equator  increases  with  Reynolds  numbers,  the  thickness  of  the  region  where 
a  general  outward  flow  from  the  equator  can  be  discussed  decreases. 


0.12 


u  0.04  - 


Figure  7.  Variation  of  Radial  Velocity  with  0  for  Re  =  200 
-.-.,  r=  1 .0543;  — ,  r=  1 . 1 1 29; ,  r=  1 .  165 1 

Figure  8  shows  the  variation  of  radial  velocity  with  radial  distance  at  the  equator 
for  different  values  of  Reynolds  numbers.  Details  of  this  swirling  jet  at  the  equator  were 
given  in  [Ref.  27]  and  [Ref.  28]. 
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Figure  8.  Variation  of  Radial  Velocity  with  Radial  Distance  at  0=7t/2  o,  Re=20 

+,  Re=200;  *,  Re=3000 

Figure  9  shows  the  variation  of  dimensionless  vorticity  at  the  surface  with  0  at 
different  Reynolds  numbers.  The  interest  of  this  figure  is  to  show  that  no  separation 
occurs  over  the  surface  of  the  sphere,  since  a  sign  change  must  take  place  as  a  point  of 
separation  is  passed. 


10  20  30  40  50  60  70  80  90 


Figure  9.  Variation  of  Dimensionless  Vorticity  with  6 
o,  Re=20;  +,  Re=200;  *,  Re=3000 
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It  is  known  that  the  boundary  layer  solution  to  be  approached  as  Reynolds 
numbers  increases,  especially  near  the  poles,  and  it  is  accurate  to  0(Re'1/2).  Therefore,  we 


might  expect  that  Re 


_1/2  do> 


at  r  =  1  and  0  =  0  will  be  of  the  form 


d(x) 


Re,"Haejs:. 


nhS-L  =A  +  BRe, 


-1/2 


(5.7) 


where  A  and  B  are  constants. 


1/2  d(0 


If  Re    ~  -r—  at  r  =  1  and  0  =  tt/2  is  evaluted,  it  can  be  seen  that  this  function  is  an 
ou 


increasing  function  of  Re.  It  can  be  assumed  that  this  expression  is  in  the  form  of 


1 


Re 


1/2 


'da) 


=  CReaa 

r=l  a 

6=71/2 


(5.8) 


where  C  and  a  are  constants. 

The  values  of  A,  B,  C  and  a  were  found  by  Dennis,  Ingham  and  Singh  [Ref.  26]  to 
be  0.51,  0.56,  0.6  and  V*  respectively.  The  modified  values  of  B  and  C  for  Reynolds 
numbers  based  on  diameter  are  0.79  and  0.5  respectively. 


Red 

1/5    ^03 
oo    e=o 

0.79 

0.51  + rpr 

Re/2 

,/5  dco 

Red-"-fc-)r=1 

0.5Red1/4 

10 

0.695 

0.760 

0.43 

0.89 

20 

0.639 

0.687 

0.806 

1.06 

50 

0.597 

0.622 

1.15 

1.33 

100 

0.576 

0.589 

1.50 

1.58 

200 

0.552 

0.566 

1.79 

1.88 

1000 

0.528 

0.535 

2.81 

2000 

0.523 

0.528 

3.34 

Table  3.  Calculated  Values  of  Re""2  —  at  r  =  1 ,  0  =  0  and  r  =  1 ,  0  =  n/2 

ou 
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VI.  CONCLUSION  AND  RECOMMENDATIONS 


In  addition  for  the  limiting  cases  of  Re«l  and  Re»l  for  which  Stokes  flow  and 
boundary  layer  solutions  are  available,  the  numerical  results  showed  very  good  agreement 
with  the  those  previously  obtained  by  other  numerical  techniques.  By  changing  some  of 
the  parameters  such  as  n*  and  L,  the  agreement  with  previous  numerical  results  could  be 
refined,  although  there  is  no  proof  that  the  other  studies  are  correct  within  three  decimal 
digits. 

Choosing  the  mapping  parameter,  L,  should  be  done  carefully  so  as  to  have  many 
grid  points  in  the  regions  of  large  gradients.  A  much  more  powerful  mapping  should  be 
determined  that  provides  a  finer  grid  in  two  completely  different  regions,  namely  the 
boundary  layer  region  and  the  re-circulation  region  (which  though  spurious,  must  be 
captured  accurately  by  the  numerical  method  for  other  calculations). 

A  different  FFT  technique  can  be  used  that  allows  a  number  of  grid  points  in  the  0 
direction  other  than  the  ones  that  are  a  power  of  two.  This  would  provide  a  relaxation  in 
the  choice  of  the  number  of  grid  points  for  different  values  of  the  Reynolds  number  which 
greatly  affect  the  memory  requirement.  The  problem  is  solved  for  0  =  0  to  9  =  n  with  the 
aid  of  sine  expansions.  Since  the  problem  is  also  symmetric  with  respect  to  7t/2,  by  taking 
care  of  the  even  function  properties  it  can  only  be  solved  for  one  quarter  which  simply 
doubles  the  grid  points  in  0  direction  with  the  same  memory  requirement. 

The  Green's  Function  Method  that  deals  with  the  Neumann  boundary  conditions 
should  not  be  considered  a  trivial  detail  and  must  be  used  to  have  a  stable  convergence 
process. 

Obtained  modified  stream  function  values,  C,  can  be  easily  used  to  solve  the  heat 
transfer  problem  due  to  convection  from  a  spinning  sphere.  Derived  additional  energy 
equations  are  as  follow 
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T(t  +  At) 


~ 

~>  ~ 

2 

r2D2 

r 
-Ec— 
At_ 

=  -T(t) 

r2D2  +Ec— 

At 

+  Ec~[3K(t)-K(t-At)]       (6.1) 


where  Ec  =  Pr  Re. 


The  nonlinear  term,  K,  in  Equation  (6.1)  consists  of  5  parts 


dTdC 
36  dr 


dTdC 
dr  96 


CrCot{6)  — 
dr 


-2Ec 


Sin2(d) 


Since,  for  forced  convection  only,  Equation  (6. 1)  is  uncoupled  from  Navier  Stokes 
equations,  the  steady  state  values  of  modified  stream  function,  C,  can  be  used  in  the 
energy  equation  after  solving  the  fluid  problem  which  decreases  the  additional  memory 
requirement  due  to  different  block  matrixes  in  Equation  (6. 1). 
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APPENDIX  A.  PROGRAM  STRUCTURE 


In  the  first  part  of  the  program,  problem  parameters  were  defined  such  as  Re,  dt, 
L,  r_inf,  M,  N.  Also  calculations  that  are  not  dependent  on  time  were  performed  in  this 
part.  To  gain  some  additional  memory  Dl,  D2  and  D3  were  calculated  with  D_Bmat 
function  and  saved  as  diagonal  elements  only.  Non-diagonal  elements  of  these  matrixes 
were  included  in  each  iteration. 

Calculations  showed  that  the  matrix  Ml  (derivative  operator  in  r  direction)  and 
consequently  M3  and  block  matrixes  that  are  diagonals  of  Dl,  D2  and  D3  are  ill- 
conditioned.  To  avoid  numerical  errors  in  solving  equation  systems,  the  Singular  Value 
Decomposition  method  was  used.  Necessary  SVD  decomposition  of  the  matrixes  and 
calculations  of  co  and  c  that  are  used  in  the  Green's  Function  Method  were  performed  in 
this  part.  The  built-in  function  in  MATLAB,  svd,  was  used  for  decomposition. 

Time  dependent  calculations  basically  consist  of  four  steps.  In  the  first  step,  u$ 
was  calculated  from  Equation  (3.32).  In  the  second  step,  particular  solutions  of  co  and  c 
were  found  from  Equations  (3.25)  and  (3.26).  In  the  third  step,  H,  which  is  used  in  the 
Green's  Function  Method,  was  found  and  finally  homogenous  and  the  particular  solutions 
were  combined  to  get  co  and  c. 

Most  of  the  calculations  were  performed  in  spectral  domain  which  is  a  feature  of 
the  numerical  technique.  Only  nonlinear  operations  were  done  in  physical  domain. 
Nonlinear  terms  were  calculated  with  the  functions  nonlin  and  nonlinsp  for  co-c  and  u 
equations  respectively. 

A  convergence  criterion  was  not  used  in  the  program.  For  high  Reynolds  Number 
(above  200),  some  over-shooting  occurs  before  steady  state.  This  situation  might  have 
caused  wrong  results  if  a  criterion  was  used.  Instead,  the  results  were  checked  during 
calculations  and  the  program  was  terminated  after  seeing  steady  state. 
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Solsvd  :  This  function  was  used  to  solve  equations  after  SVD  decomposition.  The 
limit  used  in  this  function  reduces  numerical  errors  due  to  ill-condition  matrixes  by 
skipping  close  to  singular  values  in  each  iteration. 

Solsvd  1  :  Since  this  function  was  used  for  solving  homogenous  parts  in  the 
Green's  Function  Method  and  performed  once  before  the  time-loop,  no  limit  was  used  in 
solsvd  1. 

Phy  :  This  function  was  used  for  transformation  of  the  variables  from  spectral 
domain  to  physical  domain 

Operat  :  Derivatives  in  0  direction  in  the  spectral  domain  were  performed  with  the 
function  operat.  Since  these  terms  were  used  only  in  nonlinear  terms,  the  output  from  this 
function  is  in  the  physical  domain. 

Variables  that  were  used  in  the  program  are 

L         :  Mapping  parameter. 

Re       :  Reynolds  Number. 

dt        :  Time  interval  between  iterations. 

r_inf    :  Outer  boundary. 

step     :  Number  of  iterations. 

N        :  Number  of  grid  points  in  0  direction. 

M        :  Number  of  grid  points  in  r  direction. 

xcor,  ycor  :  Cartesian  coordinates  of  grid  points  with  respect  to  the  center  of  the 

sphere  (used  for  the  presentation  of  the  results). 

Ml      :  Derivative  operator  in  r  direction. 

Dl,  D2,  D3  :  Matrixes  defined  as  in  Equation  (4.1 1). 

w         :  Vorticity. 

c  :  Potential  function  as  in  Equation  (3.15). 

u  :  Velocity  component  in  <I>  direction. 

wtil,  ctil  :  Homogenous  part  of  the  a)  and  c  defined  as  in  the  Green's  Function 

Method. 

Aij,  H,  landa  :  Variables  defined  as  in  the  Green's  Function  Method. 
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wp,  cp  :  Particular  solution  of  co  and  c  defined  as  in  the  Green's  Function  Method, 
torq      :  Dimensionless  torque  defined  as  in  Equation  (5.2). 
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APPENDIX  B.  PROGRAM  CODES 


%%  define  program  parameters 

clear  all 

global  EMI  M2M3 

L=7;Re=2;dt=0.05; 

r_inf=150;step=4500; 

%%  define  grid  points 

N=32;n=[0:N-l]'; 

teta=n.*pi./N; 

M=32;m=[0:M]'; 

z=cos(pi.*m./M); 

b=l+2*L/(r_inf-l); 

r=l+L.*(l+z)./(b-z); 

global  N  M 

xcor=cos(teta)*r'; 

ycor=sin(teta)*r'; 

E=emat(M); 

disp('E  matrix  is  done') 

%%  define  Ml, M2,M3 

dz_dr=(b.*(L+r- 1  )-b.*(r- 1  )+L)./(L+r- 1  ).A2; 

B=diag(dz_dr); 

M1=B*E; 

R=diag(r);global  R  Re  dt 

M2=R*M1; 

M3=R*(2.*eye(M+ 1  )+M2)*M  1 ; 

%%  define  Dl  D2  D3 

[D 1  ,D2,D3]=D_Bmat(M,N); 

%  define  initial  w  c  wp  cp  u 

wp=zeros(M+ 1  ,N);cp=wp;  w=wp;c=cp; 

dum  1  =sin(teta)*(  1 7r.A2)';dum2=[dum  1  ;zeros(  1  ,M+ 1 );- 1  .*flipud(dum  1  (2:N,: ))] ; 

dum3=fft(dum2); 
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u=wp; 


%  define  kmat  (for  teta  der.)  and  cotteta  (for  nonlinear  terms) 

dum=ones(l,M+l); 
kmat=n*dum; 

cotteta=[0;cot(teta(2:N))]*dum; 
global  kmat  cotteta  teta 

%  svd  decomposition  for  Dl  D3  and  define  rsquare  (for  nonlinear  terms) 

Dlu=zeros(N*(M+l),M+l); 

Dll=Dlu;Dlp=Dlu;D31=Dlu;D3u=Dlu;D3p=Dlu; 

rsquare=zeros(M+ 1  ,N); 

fori=l:N 
bas=(i-l)*(M+l)+l; 
son=bas+M; 
rsquare(:,i)=r.A2; 

Dl(bas,:)=[l  zeros(l,M)];  %  modify  Dl  for  w  @  r=inf 
Dl(son,:)=[zeros(l,M)  1];  %  modify  Dl  for  w  @  r=l 
[D 1  l(bas:son,:),D  1  u(bas:son,:),D  1  p(bas:son,:)]=svd(D  1  (bas:son,:)); 
D3(bas,:)=[l  zeros(l,M)];  %  modify  D3  for  c  @  r=inf 
D3(son,:)=[zeros(l,M)  1];  %  modify  D3  fore  @  r=l 
[D31(bas:son,:),D3u(bas:son,:),D3p(bas:son,:)]=svd(D3(bas:son,:)); 

end 

global  rsquare 

%  define  c_  c_teta  c_r  u_  u_teta  u_r 

c_=zeros(N,M+ 1  );c_teta=c_;c_r=c_; 
u_=l.*sin(teta)*(l./r.A2)'; 
u_teta=  1  .*cos(teta)*(  1  ./r.A2)'; 
u_r=-2.*sin(teta)*(  1  ./r.A3)'; 
global  c_  c_teta  c_r  u_  u_teta  u_r 

%  define  c  boundary  at  r=l 

cbrl=zeros(l,N); 

%  begin  green  function  wtil  ctil 

wtil=zeros(M+ 1  ,NA2);ctil=wtil; 
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forjj=l:N 
RR=zeros(M+l,N); 
duml=[zeros(l,jj-l)  1  zeros(l,N-jj)]; 
dum2=[duml  0  -l.*fliplr(duml(2:N))]; 
dum3=fft(dum2); 
dum4=imag(dum3(  1  :N)); 
RR(M+l,:)=dum4; 
wtil(:,jj*N)=solsvd  1  (D 1 1((N- 1  )*(M+ 1 )+ 1  :(N- 1  )*(M+ 1 )+ 1 +M,:),. 

Dlu((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

Dlp((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

RR(:,N)); 

fori=N-l:-l:l 
sum=zeros(M+ 1,1); 
ii=(jj-l)*N+i; 
forj=i+l:N 

if  rem(i+j,2)==2 
sum=sum+(-2*(i-l)).*wtil(:,(jj-l)*N+j); 

end 
sum(l)=0; 
sum(M+l)=0; 
end 
RR_=RR(:,i)-sum; 

wtil(:,ii)=solsvdl(Dll((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlp((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
RRJ; 

end 

tempwtil=(-l.*rsquare).*wtil(:,(jj-l)*N+l:jj*N); 
tempwtil(M+l  ,:)=zeros(  1  ,N); 
tempwtil(  1  ,:)=zeros(  1  ,N); 
ctil(:jj*N)=solsvdl(D31((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),.. 

D3u((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

D3p((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

tempwtil(:,N)); 

fori=N-l:-l:l 
ii=(jj-l)*N+i; 
sum=zeros(M+ 1,1); 
forj=i+l:N 

if  rem(i+j,2)==2 
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sum=sum+(-2*(i-l)).*ctil(:,(jj-l)*N+j); 
end 
sum(l)=0; 
sum(M+l)=0; 
end 
wtil_=tempwtil(:,i)-sum; 

ctil(:,ii)=solsvdl(D31((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),. 
D3u((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
D3p((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
wtilj; 
end 
end 

%  find  Aij  matrix 

forjj=l:N 

blok=ctil(:,(jj-l)*N+l:jj*N); 

duml=phy(blok')'; 

dum2=(Ml*duml)'; 

dum3=[dum2;  zeros(l,M+l) ;  -l.*flipud(dum2(2:N,:))]; 

dum4=fft(dum3);dum5=imag(dum4(  1  :N,:));dum5=dum5'; 

ctilr_l(jj,l:N)=dum5(M+l,:); 
end 

Aij=ctilr_l'; 

%  find  svd  of  Aij  for  landa  solving 

Aij(l,:)=[l  zeros(l,N-l)]; 
[Aiju,Aijs,Aijv]=svd(Aij); 

%  begin  time  loop 

for  ops=l:step 

if  ops==l  ;prc=c;prw=w;pru=u;end 

%  findRl 

Rl=zeros(M+l,N); 

Rl(:,N)=D2((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:)*w(:,N); 
for  i=N-l:-l:l 
summ=zeros(M+ 1,1); 
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for  j=i+l:N 

if  rem(i+j,2)==2 

summ=summ+(-2*(i-l)).*w(:,j); 
end 
end 

Rl(:,i)=D2((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)*w(:,i)+summ; 
end 
R1=-1.*R1; 

R 1  sp=zeros(M+ 1  ,N); 

Rlsp(:,N)=D2((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:)*u(:.N); 
for  i=N-l:-l:l 
summ=zeros(M+ 1,1); 
forj=i+l:N 

if  rem(i+j,2)==2 

summ=summ+(-2*(i- 1  )).*u(:,j); 
end 
end 

Rlsp(:,i)=D2((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)*u(:,i)+summ; 
end 
Rlsp=-l.*Rlsp; 

%    findR2sp 

R2asp=nonlinsp(u,c); 
R2bsp=nonlinsp(pru,prc); 
R2sp=3/2.*R2asp-l/2.*R2bsp; 
R2sp=(- 1  *Re).*rsquare.*R2sp; 

Rtotsp=Rlsp+R2sp; 

%  find  new  u 

%%  apply  BC  for  u 

Rtotsp(  1  ,:)=zeros(  1  ,N);    %  BC  for  u=0  at  r_inf 
Rtotsp(M+l  ,:)=zeros(  1  ,N);    %  BC  for  u  at  r=l 

u(:,N)=solsvd(Dll((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
Dlu((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
Dlp((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
Rtotsp(:,N)); 

for  i=N-l:-l:l 
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sum=zeros(M+ 1,1); 
forj=i+l:N 

if  rem(i+j,2)==2 
sum=sum+(-2*(i-l)).*u(:,j); 

end 
sum(l)=0; 
sum(M+l)=0; 
end 
Rtotsp_=Rtotsp(:,i)-sum; 

u(:,i)=solsvd(Dll((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 

Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlp((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Rtotsp_); 

end 

%    findR2 

R2a=nonlin(w,c,u); 

R2b=nonlin(prw,prc,pru); 

R2=3/2.*R2a-l/2.*R2b; 

R2=(- 1  *Re).*rsquare.*R2; 

Rtot=Rl+R2; 

prw=w;prc=c;pru=u; 

%%  green  function  solution 

%  find  new  wp 

%%  apply  BC  for  wp 

Rtot(  1  ,:)=zeros(  1  ,N);    %  BC  for  wp=0  at  r_inf 
Rtot(M+l,:)=zeros(l,N);  %  BC  for  wp=0  at  r_l 

wp(:,N)=solsvd(D  1 1((N- 1  )*(M+ 1 )+ 1  :(N- 1  )*(M+ 1 )+ 1 +M,:),. 
Dlu((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
Dlp((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
Rtot(:,N)); 

for  i=N-l:-l:l 
sum=zeros(M+ 1,1); 
forj=i+l:N 

if  rem(i+j,2)==2 
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sum=sum+(-2*(i- 1  )).*  wp(:,j); 
end 
sum(l)=0; 
sum(M+l)=0; 
end 
Rtot_=Rtot(:,i)-sum; 

wp(:,i)=solsvd(Dll((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlp((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
RtotJ; 

end 

%  solve  for  new  cp 

tempwp=-l  .*rsquare.*wp; 

tempwp(M+l,:)=cbrl ;  %  BC  for  cp=-c_r  at  r_l  zero  for  pure  spin 

tempwp(l,:)=zeros(l,N);  %  BC  for  cp=0  at  r_inf 

cp(:,N)=solsvd(D31((N- 1  )*(M+ 1 )+ 1  :(N- 1  )*(M+ 1 )+ 1 +M,:),... 
D3u((N- 1  )*(M+ 1 )+ 1  :(N- 1  )*(M+ 1 )+ 1 +M,:),... 
D3p((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
tempwp(:,N)); 

for  i=N-l:-l:l 

sum=zeros(M+ 1,1); 
forj=i+l:N 

if  rem(i+j,2)==2 
sum=sum+(-2*(i-l  )).*cp(:,j); 

end 
sum(  1  )=0; 
sum(M+l)=0; 
end 
w_=temp  wp( : ,  i)-sum; 

cp(:,i)=solsvd(D31((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
D3u((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
D3p((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
w_); 
end 

%  find  landa 
%  find  H 
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duml=phy(cp')'; 

dum2=(Ml*duml)';  %  NxM+1 

dum3=- 1  .*c_r-dum2; 

dum4=[dum3;zeros(l,M+l);-l.*flipud(dum3(2:N,:))]; 

dum5=fft(dum4); 

dum6=imag(dum5(  1  :N,:)); 

dum6=dum6';  %  M+lxN 

H=dum6(M+l,:);H=H';  %  Nxl 

H(1)=0; 

landa=solsvd(Aiju,Aijs,Aijv,H); 

%  find  complete  w  &  c 

dum  1  c=zeros(M+ 1  ,N);dum  1  w=dum  1  c ; 
forjj=l:N 

blokc=ctil(:,(jj-l)*N+l:jj*N); 

blokw=wtil(:,(jj-l)*N+l:jj*N); 

dum2c=landa(jj).*blokc; 

dum2w=landa(jj).*blokw; 

dum  1  c=dum  1  c+dum2c; 

dum  1  w=dum  1  w+dum2w; 
end 

w=wp+dum  1  w; 
c=cp+dumlc; 

dum  1  =(phy  (u')+u_)'; 

dum2=M  1  *dum  1  -dum  1 ; 

dum3=dum2(M+l,:)'.*(sin(teta).A2); 

torq(ops)=-8*pi/Re*trapz(teta,dum3); 

if  ops>l 

err(ops)=torq(ops)-torq(ops- 1 ); 

end 

save  tez4.mat  ucwRe  ops  r_inf  torq  L  N  M  dt; 

end  %  end  of  time  loop 

function  [Dl,D2,D3]=D_Bmat(M,N) 

global  M3  R  dt  Re 

k=(M+l)*N; 

Dl=zeros(k,M+l); 
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D2=zeros(k,M+l); 
D3=zeros(k,M+l); 
fori=l:N 
forj=i:N 
ifi=j 
Dl((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 
=M3+(- 1  *(i- 1  )*(i).*eye(M+ 1  )-Re/dt.*R  A2); 

D2((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 
=M3+(-l*(i-l)*(i).*eye(M+l)+Re/dt.*R.A2); 

D3((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 
=M3+(-l*(i-l)*(i).*eye(M+l)); 

end 
end 
end 

function  E=emat(M) 
E=zeros(M,M); 
m=0:M; 

z=cos(pi.*m./M); 
for  i=0:M 
forj=0:M 
if  i~=j 
if  ((i==0  I  i=M)  &  (j==0  I  j==M))  I  ((i~=0  &  i~=M)  &  (j~=0  &  j~=M)) 

c=l;else 
if  (i==o  I  i==M)  &  (j~=0  &  j~=M) 

c=2;else 
if  (i~=0  &  i~=M)  &  (j==0  I  j==M) 
c=l/2;end;end;end; 

E(i+lj+l)=c*(-l)A(i+j)/(z(i+l)-z(j+D); 
else; 
if  i==0 

E(i+l,j+l)=(2*MA2+l)/6; 
else 
if  i==M 

E(i+l,j+l)=(2*MA2+l)/(-6); 
else 

E(i+ 1  ,j+ 1  )=- 1  *z(j+ 1  )/(2*(  1  -z(j+ 1  )A2)); 
end 
end 
end 
end 
end 
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function  R2=nonlin(w,c,u) 

global  kmat  rsquare  cotteta  N  M  Ml  c_  c_teta  c_r  teta  u_  u_r  u_teta 

%  variable  list 

% 

%  dw/dteta=non  1  dc/dteta=non2  dw/dr=non3  dc/dr=non4  (in  phy) 

%  du/dteta=non5  du/dr=non6 

%  transpose  matrixes  for  colomn  operations  in  teta  direction 

r=sqrt(rsquare)';w=w';c=c';u=u'; 

nonl=operat(w); 

non2=operat(c); 

non5=operat(u); 

wph=phy(w);cph=phy(c);uph=phy(u); 

non3=[Ml*wph']'; 

non4=[Ml*cph']';non4(:,M+l)=-0.5.*sin(teta); 

non6=[Ml*uph']'; 

terml=-l./r.*non3.*(non2+c_teta); 

term2=  1  ./r.*non  1  .*(non4+c_r); 

term3=- 1  .*cotteta./r.*wph.*(non4+c_r); 

term4=-l.*cotteta./r.*non3.*(c_+cph); 

term5=wph./rsquare'.*(non2+c_teta); 

term6=(c_+cph)./rsquare'.*non  1 ; 

term7=l.*cotteta.*(non6+u_r).*(uph+u_).*2./r; 

term8=-l.*(non5+u_teta).*2.*(uph+u_)./rsquare'; 

termO=term  1  +tenri2+teiTn3+term4+term5+term6+term7+terrn8 ; 

ter=[termO;zeros(l,M+l);-l.*flipud(termO(2:N,:))]; 

RR2=fft(ter); 

R2=imag(RR2(  1  :N,:));R2=R2'; 

function  R2sp=nonlinsp(u,c) 

global  kmat  rsquare  cotteta  N  M  M 1  c_  c_teta  c_r  teta  u_  u_r  u_teta 

%  variable  list 

% 

%  du/dteta=non  1  dc/dteta=non2  du/dr=non3  dc/dr=non4  (in  phy) 

%  transpose  matrixes  for  column  operations  in  teta  direction 

r=sqrt(rsquare)';u=u';c=c'; 

nonl=operat(u); 

non2=operat(c); 

uph=phy(u);cph=phy(c); 

non3=[Ml*uph']'; 

non4=[Ml*cph']'; 

terml=-l./r.*(non3+u_r).*(non2+c_teta); 

term2=  1 7r.*(non  1  +u_teta).*(non4+c_r); 

term3=cotteta./r.*(uph+u_).*(non4+c_r); 

term4=-l.*cotteta./r.*(non3+u_r).:ic(c_+cph); 
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term5=-l.*(uph+u_)./rsquare'.*(non2+c_teta); 

term6=(c_+cph)./rsquare'.*(nonl+u_teta); 

termO=term  1  +teim2+teim3+term4+term5+teim6; 

ter=[termO;zeros(  1  ,M+ 1 );- 1  .*flipud(termO(2:N,:))]; 

RR2=fft(ter); 

R2=imag(RR2(  1  :N,:));R2sp=R2'; 

function  y=operat(x) 

global  N  M  kmat 

re=zeros(N,M+l); 

X=re+i.*x; 

Y=i.*kmat.*X; 

Y_ful=[Y;zeros(l,M+l);flipud(Y(2:N,:))]; 

y=ifft(Y_ful); 

y=real(y(l:N,:)); 

function  y=phy(x) 

global  N  M 

re=zeros(size(x)); 

X=re+i.*x; 

Y_ful=[X;zeros(  1  ,M+ 1  );conj(flipud(X(2:N,:)))]; 

dum=ifft(Y_ful); 

y=real(dum(l:N,:)); 

function  y=solsvd(u,s,v,b) 
w=diag(s); 

wmin=max(  w)*  1  e- 1 2; 
for  i=l:length(w) 

if  w(i)<wmin;ww(i)=0;else;ww(i)=l/w(i);end 
end 

www=diag(ww); 
y=v*www*(u'*b); 

function  y=solsvdl(u,s,v,b) 

w=diag(s); 

for  i=l:length(w) 

ww(i)=l/w(i); 
end 

www=diag(ww); 
y=v*www*(u'*b); 
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